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A model for two-dimensional colloids confined laterally by "structured boundaries" (i.e., ones that 
impose a periodicity along the slit) is studied by Monte Carlo simulations. When the distance D 
between the confining walls is reduced at constant particle number from an initial value Do, for which 
a crystalline structure commensurate with the imposed periodicity fits, to smaller values, a succession 
of phase transitions to imperfectly ordered structures occur. These structures have a reduced number 
of rows parallel to the boundaries (from nton— Iton — 2 etc.) and are accompanied by an almost 
periodic strain pattern, due to "soliton staircases" along the boundaries. Since standard simulation 
studies of such transitions are hampered by huge hysteresis effects, we apply the phase switch Monte 
Carlo method to estimate the free energy difference between the structures as a function of the misfit 
between D and Do, thereby locating where the transitions occurs in equilibrium. For comparison, 
we also obtain this free energy difference from a thermodynamic integration method: the results 
agree, but the effort required to obtain the same accuracy as provided by phase switch Monte Carlo 
would be at least three orders of magnitude larger. We also show for a situation where several 
"candidate structures" exist for a phase, that phase switch Monte Carlo can clearly distinguish the 
metastable structures from the stable one. Finally, applying the method in the conjugate statistical 
ensemble (where the normal pressure conjugate to D is taken as an independent control variable) 
we show that the standard equivalence between the conjugate ensembles of statistical mechanics is 
violated. 



I. INTRODUCTION 



Periodically ordered arrays of nanoparticles, colloidal 
crystals, crystalline mesophases formed from surfactant 
molecules or block copolymers, etc. are all examples of 
complex periodic structures that can occur in soft mat- 
ter systems. Since often the interactions between the 
constituent particles of these structures are to a large 
degree tunable, one has the possibility of producing ma- 
terials vifith "tailored" properties vifhich have potential 
applications in nanotechnological devices ^^. When 
seeking to provide theoretical guidance for understanding 
structure-property relations in such complex soft mat- 
ter systems, a basic issue is how to judge the relative 
stability of competing candidate structures, i.e. to dis- 
tinguish the stable structure (having the lowest free en- 
ergy) from the metastable ones. For standard crystals 
formed from atoms or small molecules, this question can 
be answered by comparing ground state energies of the 
competing structures (and -if necessary- also taking en- 
tropic contributions from lattice vibrations into account, 
within the framework of the harmonic approximation). 
In soft matter systems, disorder in the structure and ther- 
mally driven entropic effects rule out such an approach, 
and hence there is a need for computer simulation meth- 
ods that compute the free energy of the various complex 
structures. However, as is well known, the free energy of 
a model system is not a direct output of either Molec- 
ular Dynamics or Monte Carlo simulations, and special 
techniques have to be used [g-llH. 

In principle, one can obtain the absolute free energy 
of a structure by linking it to some reference state of 



known free energy by means of thermodynamic integra- 
tion (TI) Wldl- The strengths of TI are that it is both 
conceptually simple and often straightforward to imple- 
ment. Its principal drawback is that the quantity of in- 
terest, namely the free energy difference between candi- 
date structures is typically orders of magnitude smaller 
than the absolute free energies of the individual struc- 
tures which TI measures. Essentially, therefore, TI es- 
timates a small number by taking the difference of two 
large ones; As a consequence, the precision of the method 
is limited and an enormous (even sometimes wasteful) in- 
vestment of computer resources rnay be needed to resolve 
free energy difference accurately [3|. 

A much more elegant approach, albeit one which is not 
quite so easy to implement as TI, is the "phase switch 
Monte Carlo" [17H23| technique. This method is poten- 
tially more powerful than TI because it focuses directly 
on the small free energy difference between the structures 
to be compared, rather than their absolute free energies. 
In previous work, the precision of the method was demon- 
strated in the context of measurements of the free energy 
difference between fee and hep structures of hard spheres 
|20| , the phase behaviour of Lennard- Jones crystals [20| 
and as a means of studying liquid-solid phase transitions 
[l8{ . In the latter case, simple model systems contain- 
ing only a few hundred particles could be studied, while 
for the study of the fcc-hcp free energy difference [13, l2l| 
larger systems of up to 1728 particles could be studied by 
virtue of the fact that these crystals only differ in their 
packing sequence of close-packed triangular defect-free 
lattice planes. However, it is an open question as to what 
system sizes one can attain with the phase switch method 
for more general crystalline systems, including - as in the 



present work ~ ones which exhibit considerable structural 
disorder ("soliton staircases", see below). Furthermore, 
there have hitherto been no like-for-like comparisons of 
the TI and phase switch methods on the same system, so 
whilst their are good reasons for presuming the superior- 
ity of phase switch (in terms of precision delivered for a 
given computational investment), this has never actually 
been quantified. 



In the present paper, we address these matters, con- 
sidering as a generic example a two-dimensional colloidal 
crystal in varying geometrical confinement [2J-|28| . As is 
well-known, two dimensional colloidal crystals are exper- 
imentally much studied model systems 29-40] compris- 
ing, for example, polystyrene spheres containing a super- 
paramagnetic core adsorbed at the air- water-interface. 
Applying a magnetic field oriented perpendicular to this 
interface creates a repulsive interaction that scales like 
r~^, (r being the particle separation), whose magnitude 
is controlled by the magnetic field strength [29|. Lat- 
eral confinement of such two-dimensional crystals can 
be effected mechanically or by laser fields (if the latter 
are also applied in the bulk of such a crystal, one can 
study laser- induced melting and/or freezing [4l| - |44| ). Of 
course, there exist many related problems in rather dif- 
ferent physical contexts ("dusty plasmas" pa . |46| . i.e. 
negatively charged Si02 fine particles with 10/zm diam- 
eter in weakly ionized rf discharges; lattices of confined 
spherical block copolymer micelles [43]; vortex matter 
in slit channels i48|, etc.). However, our study does not 
address a specific system, rather we focus on the method- 
ological aspects of how one can study such problems by 
computer simulation. 



The outline of the present paper is a follows. In Sec. 
2, we summarize the key facts about our model, namely 
strips of two-dimensional crystals confined between two 
walls where structural phase transitions may occur when 
the distance between the (corrugated) rigid boundaries 
is varied |25l - [28l |49| - |5]| (i.e., a succession of transitions 
in the number of crystal rows n parallel to the walls oc- 
cur, n— )-n — 1— >n— 2, with increasing compression, 
accompanied by the formation of a "soliton staircase" 
at the walls |25l428| ). In Sec. 3, the methods that are 
used are briefly described: the thermodynamic integra- 
tion method of Schmid and Schilling [1^, |lg] is used as a 
baseline, while the main emphasis is on the phase switch 
Monte Carlo method (implementation details of which 
are deferred to an Appendix). In Sec. 4 we describe 
the results of the application of these techniques to the 
model of Sec. 2. We show that phase switch Monte Carlo 
p^l - l2g ] can accurately locate the phase transitions de- 
spite the need to deal with thousands of particles, and is 
orders of magnitude more efficient than thermodynamic 
integration. Sec. 5 summarizes some conclusions. 



II. STRUCTURAL TRANSITIONS IN 

CRYSTALLINE STRIPS CONFINED BY 

CORRUGATED BOUNDARIES: 

PHENOMENOLOGY 

Here we introduce the model for which our methodol- 
ogy is exemplified, and recall briefly the main findings 
concerning the rather unconventional transitions that 
have been detected |25l - l28| . as far as they are relevant 
for the present study. 

We consider monodisperse colloidal particles in a 
strictly two-dimensional geometry, which then are 
treated like point particles in a plane interacting with 
a suitable effective potential V{r) that depends only on 
the interparticle distance r. In the real systems [23, l3ll - 
1351 ] this potential is purely repulsive, but due to the 
magnetostatic dipole-dipole interaction (whose strength 
is controlled by the external magnetic field) it is very 
slowly decaying, V{r) ex r~^. Since we here are not con- 
cerned with quantitative comparisons with real experi- 
mental data on such systems, we simplify the problem 
by adopting a computationally more efficient r~^^ po- 



tential, in accord with previous work |25l - l28| . Moreover, 
to render it strictly short-ranged, we introduce a cutoff 
Tc, such that V(r > r^) = 0, and employ a smoothing 
function to make V{r) differentiable at r = r^. Thus, the 
model potential used is 



y(r) = e {a/rf^ - {a/r,f 



(r - rcY 



hi + ^r- r^Y 



, (1) 



with parameters Vc — 2.5a and h = O.OIct. Henceforth, 
the particle diameter a = 1 defines the length units in 
our model, and for the energy scale, e = 1 is taken, while 
Boltzmann's constant /cb = 1- It is known that at T = 
the ground state of this model is a perfect triangular 
lattice, with a lattice spacing a related to the choice of 
number density p = N/V (with N the particle number 
and V the (two-dimensional) "volume" of the system) via 



2/{V3p) 



(2) 



Assuming the physical effect of truncating the poten- 
tial can be neglected, only the choice of the combination 
X = p{e/kBTy^^ controls the phase behavior 52]. Thus, 
following previous work in the NVT-ensemblc it suffices 
to choose a single density when the temperature variation 
is considered [25, 53]. For the particular choice p = 1.05, 
the melting transition of this model is known to occur at 
about T — Tm ~ 1.35 [Sj]. Note that here we are not 
at all concerned with the peculiarities of melting in two 
dimensions [54], and hence we focus on a temperature 
deep within the crystalline phase, T — 1. Although it 
is known that the density of vacancies and interstitials 
in rf = 2 for any nonzero temperature is also nonzero 
in thermal equilibrium [54l . ISSl j. for the chosen particle 
number N = 3240 the system is essentially defect free. 
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FIG. 1. Sketch of the system geometry, showing the fixed 
wall particles (black spheres) and the mobile particles (gray 
spheres). The orientation of the coordinate axes is indicated, 
as well as the lattice spacing of the triangular lattice (a) and 
the linear dimensions Lx, D of the system. 



since the densities of these point defects at T = 1 are 
extremely small [23, Hj] . 

The geometry of the present system is a Z? x L^ 
slit, confined in the y-direction and periodic in the x- 
direction. In the y-direction there are Uy = 30 rows of 
the triangular lattice, each containing n^. = 108 parti- 
cles, stacked upon each other. The a;-direction coincides 
with a lattice direction so that L^ = Jii-a. The confin- 
ing boundaries (one at the top and one at the bottom 
of the system) each take the form of a double rows of 
particles in which the particles are rigidly fixed at the 
sites of a perfect triangular lattice (Fig. [T]). These rows 
of fixed particles represent rigid corrugated walls, essen- 
tially acting as a periodic wall potential on the mobile 
particles. While the distance of the first row at the up- 
per wall from the first row of mobile particles in the ideal 
stress-free crystal is simply D — nya^/3/2, in the follow- 
ing we are interested in the response of the system when 
the walls occur at a smaller distance, caused by a misfit 
A, defined via [Sg 



D = {ny- A)aV3/2 



(3) 



As described in the previous work [25l - l28j . standard 
Monte Carlo simulation [a, 0] allows one to study this 
model at various values of A, and also sample the stress 
"' = '^yy ~ '^xx {cap are the Cartesian components of the 
pressure tensor) applying the virial formula [a, 0] • Fig. [2] 
shows that when one starts out with the perfect crystal 
{uy = 30) with no misfit, the crystal already shows a 
small finite stress, because the rigid wall particles some- 
what hinder the vibrations of the mobile particles in their 
potential wells, but this effect is of no importance here. 
Rather we focus on the (slightly nonlinear) increase of 
the stress up to about A = Ac « 2, followed by the 
(almost) discontinuous decrease, and the subsequent in- 
creases again with further enhancement of the misfit. A 
previous structural analysis has revealed |25l - [28| that the 
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FIG. 2. Stress a plotted vs. misfit A, for a system of 
A^ — 3240 particles, and using different starting configura- 
tions having Uy = 30, Uy — 29, and Uy — 28, as indicated in 
the figure. Note the huge hysteresis of the n^ = 30 — >■ n^ = 29 
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sudden decrease of stress is due to a transition in the 
number of rows in the crystal, Uy -^ Uy — \ — 29. How- 
ever, since in the NVT ensemble the particle number is 
conserved, the n^, = 108 particles of the row that disap- 
pears must be redistributed among the remaining rows. 
A closer examination of the structure revealed that none 
of these particles enter the two rows adjacent to the rigid 
walls, instead they all go into the n^ — 3 = 27 rows of 
the system that are further away from the walls. Thus, 
in the present case, the particle number per row becomes 
n^ + n^l {ny — 3) = n^. -|- 4, and this leads to a new lattice 
spacing in the x-direction of a' = a/(l -I- A/ux), which 
is no longer commensurate with the spacing between the 
particles forming the rigid walls (or the two immediately 
adjacent layers which remain commensurate with them). 
While for the rows in the center of the system (near ny/2) 
this compression of the lattice spacing occurs uniformly 
along the a;-direction, this is not the case close to the 
walls, which provide a periodic potential (with period- 
icity a) that acts on the row of mobile particles a little 
further inside the slit. The fact that on the scale L^ 
the effective wall potential exhibits Ux minima but that 
n[^ — nx + 4: particles need to be accommodated, leads to 
the formation of a lattice of solitons close to both walls 
( "soliton staircase" ) [53, [S^ , as depicted for an idealized 
case in Fig. [3] 

In practice, the actual structure having nj, — 1 = 29 
rows that is formed in the simulations on increasing the 
misfit A beyond the critical value Ac, is generally less reg- 
ular than the 'idealized' one shown in Fig. |3l specifically, 
the relative distance between neighboring solitons showed 
a considerable variation. This comes about because (i) 
the solitons are formed from the stressed crystal with 
Uy = 30 rows via random defect nucleation events [26|, 
and (ii) the mutual interaction between neighboring soli- 
tons, which is the thermodynamic driving force towards 
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FIG. 3. a) Putting n + 1 particles in a periodic potential 
with n minima creates a soliton configuration, i.e. over a 
range of several lattice spacings particles are displaced from 
the potential minima (schematic) b) Superimposed snapshot 
pictures of 750 configurations of the particle positions, where 
for a system of Uy — 30 rows and a large misfit (A — 2.6) a 
transition to n^ — 1 = 29 rows has occured (n^; = 108 and 
T = 1.0 were chosen). The 4 solitons at each wall are visible 
due to the larger lateral displacements of the particles, leading 
to a darker region in the snapshot. Part (c) shows a close-up 
of the structure near the upper wall. Numbers shown along 
the axes indicate the Cartesian coordinates of the particles. 
Parts (b) and (c) have been adapted from Chui et al. [25|. 



a regular soliton arrangement, is very small [27[. De- 
spite this, it is nevertheless reasonable to construct "by 
hand" the expected regular structure oinx/{ny — 3) (= 4) 
solitons near each wall as a starting configuration for a 
system with 29 rows, which can subsequently be equili- 
brated [23 ■ Of course, there is no guarantee that this 
guessed structure actually is the one of lowest free en- 
ergy; but it does exhibit slightly less stress than all other 
structures that had been tested, for misfits in the range 
1.5 < A < 3, and hence has been used as a starting point 
for studies in which A was varied in this range. 

Starting from this idealized 29 row structure and de- 
creasing the misfit one finds that the 29 row structure 
is stable down to about A^ « 1.3, at which point the 
soliton lattice disappears and the system spontaneously 
transforms into a defect free structure with Uy = 30 rows 
again (Fig. [5]) . This value of A is to be compared with 
that for the reverse transition from 30 to 29 rows which 
we recall occurs at A^, « 2.0. Thus, with the standard 
Monte Carlo approach there is considerable hysteresis 
which precludes the accurate location of the transition 
point. Clearly, therefore a method is needed from which 
one can locate where the transition occurs in equilibrium. 

Similar hysteresis is observed if one starts out from the 
29 row structure but increases the misfit beyond A = 3 
(a case that has not been studied previously) . As Fig. [2] 
shows, a transition occurs to structures with Uy — 2 = 28 
rows (at about A w 4.1). Unfortunately, there seem 
to be no unique candidates for stable structures having 
ny — 2 = 28. Fig.Udisplays four candidate structures that 



FIG. 4. Configurations with A^ = 3240 particles and Uy ~2 = 
28 rows, but different configurations of the solitons. In the 
text, they are referenced as "configuration nr. 1, 2, 3, 4" from 
top to bottom. For a clear identification of the positions of 
the solitons, the method described in 231 was used. 



we have identified, each of which is at least metastable 
on simulation timescales. Depending on which of these 

28 row candidates one takes, the transition from 28 to 

29 rows on reducing the misfit occurs at anything be- 
tween A = 3.2 and 3.75. As regards the nature of the 
candidate structures, in each case 2nx = 216 extra par- 
ticles have to be distributed across the system. If we 
again keep the rows adjacent to the walls free of ex- 
tra particles, the particle number per inner row becomes 
n'^ = Ux + 2nx/{ny — A) Ki Ux + 8.3, i.e. is non-integer. If 
we kept two rows adjacent to the wall rows free of extra 
particles, we would have 9 extra particles per row, and 
thus this structure has been tried (this is configuration 
number 1 in Fig. 2]). Another structure was obtained if 
we place 4 extra particles in the rows directly adjacent 
to the walls and 8 extra particles in each of the 26 inner 
rows (configuration number 2). By energy minimization 
of a somewhat disordered structure resulting from a tran- 
sition from 29 to 28 rows a structure was obtained which 
had 9 solitons on one wall but only 8 on the other wall 
(configuration number 3). Finally another configuration 
with 8 solitons on each wall (configuration number 4) 
was found. Note that the configurations shown in Fig. |3] 
are not the actual structures at T = 1.0 but the cor- 
responding "inherent structures" found from the actual 
structures by cooling to T = 0, to clearly display where 
the solitons occur. Clearly, it again is a problem to (i) 



identiiy which of these 4 configurations with 28 rows is 
the stable one (at T = 1.0), and (ii) determine at which 
misfit the transition to the structure with 29 rows occurs. 
As we shall demonstrate below, both problems can be el- 
egantly dealt with by employing the phase switch Monte 
Carlo method. 



III. FREE ENERGY BASED SIMULATION 

METHODOLOGIES TO LOCATE TRANSITIONS 

BETWEEN IMPERFECTLY ORDERED 

CRYSTAL STRUCTURES 

A. Thermodynamic Integration 

The general strategy of thermodynamic integration is 
to consider a Hamiltonian 'H(A) that depends on a pa- 
rameter A that can be varied from a reference state (char- 
acterized by Aq) whose free energy is known, to the state 
of interest (Ai), without encountering phase transitions. 
The free energy difference AF can then be written as 



Ai 



AF = F(Ai) - F(Ao) - / dX'{d'H{\')/dX')x 



(4) 



50(a) = T2 [exp(a) - ^ e''/k\] 



k=0 



(7) 



for the choice of (j){x) written in Eq. ^. 

Then intermediate models 'H(A) to be used in Eq. (|3]) 
are chosen as 



H'(A)=Hi„t + C/rct(A) 



(8) 



where Hint describes interactions in the system, which 
then are switched on (if necessary, in several steps). The 
free energy contribution of switching on these interac- 
tions can easily be determined by a Monte Carlo simu- 
lation which includes a move that switches the interac- 
tions on and off. The logarithm of the ratio of how many 
times the states with and without interactions were vis- 
ited gives the free energy contribution. The free energy 
difference between the intermediate model where parti- 
cle interactions are turned on and potential wells are also 
turned on, and the target system with particle interac- 
tions but without potential wells, then is computed by 
thermodynamic integration, for which 



For a dense disordered system (fluid or a solid con- 
taining defects). Schilling and Schmid flU, [lal proposed 
to take as a reference state a configuration chosen at 
random from a well equilibrated simulation of the struc- 
ture of interest, at values of the external control param- 
eters for which one wishes to determine the free energy. 
Particles can be held rigidly in the reference configura- 
tion {fi '■®^} by means of a suitable external potentials. 
(We recall that a somewhat related thermodynamic in- 
tegration scheme for disordered systems known as the 
"Tethered spheres method" has already been proposed 
by Speedy [53 •) When these external potentials act, the 
internal interactions can be switched off. In practice, one 
can use the following pinning potential C/rcf(A) to create 
the reference state, where rcut is a parameter discussed 
below. 



(9H.ef (A)/9A) = (^ 0(|r, - n ■^'^V^cut)) (9) 



f/rcf(A) = A^0(|rl-' 



/'"cut) with(j){x) = X -1 



, . (5) 

Here it is to be understood that particle i is only pinned 
by well i at fl'^^, and not by other wells. However, iden- 
tity swaps need to be carried out to ensure the indis- 
tinguishability of particles. The free energy of this non- 
interacting reference system then is 



i^,et(A) = HN/V) - ln[l + {Vo/V)g4pX)] , (6) 
where /3 — [ksT]^^, Vq (in d = 2 dimensions) is Vb — 



needs to be sampled [la^ l6|. This method has been 
tested for hard spheres 1^, [3 ^ including also systems 
confined by walls from which wall excess free energies 
could be sampled |60| . 



B. Phase Switch Monte Carlo 



Trrr 



and 



The phase switch method [18l - l23l | computes directly 
the relative probabilities of two phases, by switching be- 
tween them and recording the ratio of the simulation time 
spent in each. This ratio directly yields their free energy 
difference AF via AF = ln(yl(i)/A(2)). Here A^^) and 
A^^' are the times spent in the respective phases which 
are proportional to the statistical weight of each phase 

The power of the phase switch method derives from its 
ability to leap directly from configurations of one pure 
phase to those of another pure phase (Fig. [SJ , avoiding 
the mixed phase states which - when one or both phases 
are crystalline - can be computationally problematic (see 
appendix A) . The leap is implemented as a suitable global 
Monte Carlo move. One starts out by specifying for each 
of the two phases of interest (labeled by index a = 1,2), 
a reference configuration. This can be expressed as a set 
of i = 1 . . . A^ particle positions {-Rj }• Note that the 
specific choice of a reference configuration for phase a 
does not matter (at least in principle, see Appendix), it 
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FIG. 5. Schematic comparison of (a) the standard method 
for hnking phases via a samphng path and (b) The phase 
switch method. The blobs represent the set of values of some 
macroscopic property (eg order parameter or energy) asso- 
ciated with configurations belonging to two distinct phases 
(a = 1, 2). These pure phase states (having high probability) 
are separated by a "deep valley" in the free energy landscape 
corresponding to interfacial states having a very low probabil- 
ity, (a) In the standard strategy one uses extended sampling 
to negotiate the valley, by climbing down into it from one side 
and climbing up out of it on the other, (b) The idea of phase 
switch Monte Carlo is to "jump over the valley". 



need only be a member of the set of pure phase config- 
urations that "belong" to phase a. Thus for example in 
the present case, a suitable reference configuration for the 
71 = 30 row defect-free structure could simply be a typ- 
ical configuration chosen from a simulation run on this 
structure. However, it could equally be a configuration in 
which all particles are at the lattice sites of this structure. 
Given the two reference configurations, one can express 
the position vectors r^ of each particle i in phase a as 



always be rejected by the Monte Carlo lottery. This prob- 
lem is circumvented by employing extended sampling 
methods [3, [lO, ISli that create a bias which enhances the 
occurrence of displacements {ut} for which the switch 
operation does have a sufficiently high Monte Carlo ac- 
ceptanceprobability. Such states are called "gateway 
states" |18l - [22 |: crucially, they do not need to be speci- 
fied beforehand - the system autonomously guides itself 
to them in the course of the biased sampling. 

In practice, the bias is administered with respect to 
an "order parameter" M whose instantaneous value is 
closely related to the energy cost of implementing the 
phase switch. One then introduces a weight function 
ri{M) into the sampling of the effective Hamiltonian 
which enhances the probability of the system sampling 
configurations for which the energy cost of the phase 
switch is low, thereby increasing the switch acceptance 
rate. Of course, the weight function /^(M) to be used is 
not known beforehand, and thus needs to be iteratively 
constructed in the course of the Monte Carlo sampling. 
One has a choice of ways of doing so: we have used the 
transition matrix Monte Carlo method 61 63] (see also 
the Appendix for implementation details). Alternative 
methods such as Wang-Landau sampling [64| or succes- 
sive umbrella sampling [73| could also be applied. 

Once a suitable form for the weight function /^(M) has 
been found, a long Monte Carlo run is performed, in the 
course of which both phases are visited many times. The 
statistics of the switching between phases is monitored 
by accumulating the histogram of M, which (as in all ex- 
tended sampling methods) is corrected for the imposed 
bias at the end of the simulation. Doing so yields an es- 
timate of the true equilibrium distribution P(M), which 
in general exhibits a double peaked form (one peak for 
each phase). The free energy difference between the two 
phases is simply the logarithm of the ratio of the peak 
weights as described at the start of this subsection. 

Of course, the above description was only intended to 
outline the phase switch strategy; more extensive imple- 
mentation details are given in the appendix. 



rv^"' = Rl"^ 



(10) 



where {ui} is a set of displacement vectors which mea- 
sure the deviation of each particle from the reference site 
to which it is nominally associated. Note that while there 
is a separate reference configuration for each phase, the 
single set of displacements is common to both phases. 

Let us suppose the simulation is currently in phase 
a — 1. Now the phase switch idea is to a map the current 
configuration {r^ } of this phase on to a configuration 
of phase a = 2 by switching the sets of reference sites 
from {R> '} to {R> '} but keeping the set of displace- 
ments {ui} fixed. This switch can be incorporated in a 
global Monte Carlo move. Of course, in general the set 
displacements that are typical for phase a = 1 will not be 
typical displacements for phase a — 2. As a consequence, 
in a naive implementation such a global move will almost 



IV. RESULTS 

A. Free energy differences and computational 
efficiency 

Fig. [6] shows the absolute free energies in the NVT 
ensemble for the phase with 30 rows (and no defects) 
and the phase with 29 rows and the "soliton staircases" 
(Fig. [HJj) as a function of the misfit A, as obtained from 
the thermodynamic integration method (Sec. III.l). One 
sees that these free energies are very large (note the or- 
dinate scale) and vary rather strongly with A. How- 
ever, the free energy curves with these two structures are 
barely distinct from each other, and hence a very sub- 
stantial computational effort is needed to locate, with 
meaningful accuracy, the intersection point marking the 
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FIG. 6. Absolute free energy F of systems of A'^ = 3240 par- 
ticles interacting with the potential given in Eq. ([1} in L x _D 
geometry with L — 108a, a being the lattice spacing, and 
periodic boundaries in x-direction, confined by two rows of 
fixed particles on either side in y-direction (Fig. [TJ as a func- 
tion of the misfit A Eq. (|3]). Two structures are compared: (i) 
a (compressed) triangular lattice with Uy — 30 rows contain- 
ing Ux = 108 particles per row; (ii) a lattice with Uy = 29 
rows and corresponding soliton staircases (Fig. |3]3) . 
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FIG. 7. Free energy differences between structures with 29 
and 30 rows plotted versus the misfit A. Both results obtained 
from thermodynamic integration and from the phase switch 
method are shown, as indicated. 



equilibrium transition between n = 30 and n = 29 rows. 

Fig. [7] plots the free energy difference AF versus the 
misfit, comparing the results from the thermodynamic 
integration method (points with error bars) with the re- 
sults from the phase switch method, and focusing on the 
region near the transition. One can see that within the 
errors the results of both methods agree very well with 
each other, although for the thermodynamic integration 
method the error is at least an order of magnitude larger 
than that of phase switch. We note that the predicted 
equilibrium value of the misfit at the transition point 
(At « 1.7) falls well within the hysteresis loop of Fig. [51 

Since the absolute free energies are of the order of 



20000 (for our system with N = 3240 particles) but, in 
the region of interest, free energy differences are of order 
±40 only, we have that the relative error 5F/F is of order 
1/500. Thus for thermodynamic integration, it would be 
difficult to bring the error bars down further in Fig. [T] 
The error bars for the phase switch simulation were com- 
puted from the results of four independent runs for each 
value of the misfit, and are hardly visible on the scale of 
Fig. [71 



In addition to this significant difference with respect 
to the size of the statistical errors, phase switch Monte 
Carlo also outperformed the thermodynamic integration 
method with respect to the necessary investment of com- 
puter resources. In order to obtain a suitable weight 
function for our system, at a certain value of the mis- 
fit, we let the simulation run for about 15 million steps 
(each step consisting of one sweep of local moves and 
one attempt to switch the phases). On the ZDV clus- 
ter of the University of Mainz, this takes about 4.5 days 
on a single core (though in hindsight we could have got 
away with a less smooth weight function, further reduc- 
ing the computing time of this step) . Having determined 
the weight function, we initiated four production runs 
for every value of the misfit. These runs needed again 
10 million steps each (i.e. about 3 days each) in order 
to perform a sufficient number of phase switches to yield 
results of the desired precision. Overall, then, computing 
each point of the free energy difference curve of Fig. [71 by 
phase switch took about 16.5 days of CPU time. 



In contrast to this, the thermodynamic integration 
method required a calculation not only of the free en- 
ergy difference in which we are interested, but of the free 
energy difference along the path of the thermodynamic 
integration, gradually switching off the wells of attrac- 
tion used there, and of the free energy difference between 
the state where the particle interactions were turned on 
and the state where they were turned off. This needs 
to be done for both phases separately. It is therefore 
not surprising, that considerably more CPU time was 
needed: roughly 250 days of CPU time were invested for 
each phase and for each value of the misfit to obtain the 
absolute free energy (again converting units to a single 
core). Thus, each of the 12 values of free energy differ- 
ences needed for Fig. [71 required 500 days (rather than 
16.5 days), i.e. a factor of 30 more computational effort! 
However, if we were to bring the statistical errors of the 
thermodynamic integration method a factor of 10 down 
(to make it comparable to the phase switch method), we 
would need another factor of 100 in computer time; the 
benefit of using the (clearly much more powerful) phase 
switch approach hence amounts to a gain of the order of 
lO'^ in computational resources! Of course, this is no sur- 
prise when we remember that the free energy differences 
of interest are only of the order of (1/500) of the total 
free energies for the present model system. 
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the index N will be omitted. Of course, at fixed lateral 
dimensions L a variation of D is equivalent to a variation 
of the volume V). 
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FIG. 8. Schematic description of phase transitions in thin 
films of thickness D in the conjugate NpT (left) and NVT 
(right) ensembles, for the case of a vapor to liquid transition 
(a) and the present transition where the number of rows is 
reduced (n — >■ n — 1) when either the (normal) pressure p 
increases (left) or the thickness decreases (right). Note that 
in the latter case two-phase coexistence is possible for the 
vapor-liquid transition, but not for the transition where the 
number of rows parallel to the boundaries change. For further 
explanations cf text. 



B. Ensemble inequivalence 

We turn now to a discussion of a puzzling aspect of 
the physics, namely the fact that we treat here a first- 
order structural phase transition obtained by variation of 
the distance D between the walls formed by the rigidly 
fixed particles, i.e. an extensive rather than an intensive 
thermodynamic variable. If we were concerned with the 
study of a vapor to liquid transition of a fluid in such a 
geometry, the proper way to locate a discontinuous tran- 
sition is the variation of the intensive variable thermody- 
namically conjugate to D, which is the normal pressure 
Pn (force per area acting on the walls; in the following 
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FIG. 9. a) Free energy difference AF for the transition from 
n = 30 to n = 29 rows as a function of pressure, (b) The 
distribution of the internal energy difference between the two 
phases p{E30rows — E29roms) at fixed {u}. Curves for 4 pres- 
sures near and at the transition pressure pt — 22.146 ± O.OfS 
are shown, as generated via histogram reweighting. The sim- 
ulation was run at a pressure of p = 22.13. (c) System length 
D as a function of pressure. Clearly, the curve for the stable 
phase exhibits a jump at the transition pressure. Statistical 
errors are smaller than the symbol sizes. 

To fix ideas, we remind the reader about this classical 
vapor-liquid problem in Fig. [5^): In the NpT ensemble, 



we would have a jump in volume V = LD from Vv = LDy 
(density of the vapor p„ — N/Vy) to 14 = LDg (density 
of the liquid pg = N/Vg) at the transition pressure pt- If 
we work in the conjugate NVT ensemble, of course, the 
behavior simply follows from a Legendre transform, the 
volume jump from Vv to V(, translates into a horizontal 
plateau at p = pt, and any state of this plateau is a situa- 
tion of two-phase coexistence, as schematically indicated 
in Fig. Hi). 

Of course, it is also possible to consider the present 
transition between a state of n rows to n — 1 rows in the 
NpT ensemble (Fig. [Sjj and Fig. [9];). Then it is clear 
that the transition will show up as a jump in the thick- 
ness D from Dn{= nan) to D„_i (= (n— l)a„_i), where 
a„, a„_i are the (average) distances between the lattice 
rows (or lattice planes, in three dimensional films, respec- 
tively). The corresponding phases of the n- layer state 
and (n — 1) layer state are indicated below the isotherm 
in the {p — D) plane schematically. 

However, one simply cannot construct a state of two- 
phase coexistence out of these two "pure phases" at a 
value of D intermediate between -D„-i and £>„: locally 
the n-layer state requires a thickness £>„, the (n — 1) layer 
state a thickness I?„^i, so one would have to "break" the 
walls. Of course, it is not just sufficient to have a state 
with n layers separated by a grain boundary from a state 
with (n— 1) layers at the same value of D: these domains 
are not the coexisting pure phases in the NpT ensemble! 

So the phase coexistence drawn (horizontal broken 
curve) in Fig. [8b) is unphysical, it requires a state where 
the constraining walls were broken. Requesting the in- 
tegrity of the walls is a global constraint which makes 
phase coexistence in the standard sense impossible for 
the present transitions! Thus, the rule that the differ- 
ent ensembles of statistical mechanics yield equivalent 
results in the thermodynamic limit is not true for the 
present system; in the transition region Dn-i < D < Dn 
the NVT ensemble and the NpT ensemble are not equiv- 
alent. 

Actually this is not the first time that such an ensem- 
ble inequivalence has been pointed out. A case much 
discussed in the literature is the "escape transition" of a 
single polymer chain of N beads grafted at a planar sur- 
face underneath a piston held at a distance D above the 
surface to compress the polymer 66;- 72]. For pressures 
p < Pt (where the piston is at distance -Dt,i) the chain 
is completely confined underneath the piston (which has 
the cross section of a circle in the directions parallel to 
the surface) while for p > pt the chain is (partially) es- 
caped into the region outside of where the piston acts 
(the piston distance at pT jumps to a smaller value £'t,2)- 
When we use instead D as the control variable, again a 
sharp transition occurs (for N ^ oo) at some interme- 
diate value Dt (£'t,2 < Dt < -Dt,i), since obviously it is 
simply inconceivable to have within a single chain phase 
coexistence between states "partially escaped" and "fully 
confined" , since these states are only defined via a global 
description of the whole polymer chain. 



Another case where transitions of the number n of lay- 
ers in layered structures in thin films occurs is the con- 
finement of symmetric block copolymer melts (which may 
form a lamellar mesophase of period Ap in the bulk) in 
thin films between identical walls |73l - [76| . When then 
the thickness D of such films is varied, one observes ex- 
perimentally discontinuous transitions in the number n 
of lamellae parallel to the film [74, 75] . However, when 
one considers block copolymer films on a substrate and 
does not impose the constraint of a uniform thickness but 
rather allows the upper surface to be free, then indeed 
mixed phase configurations of a region where n — 1 layers 
occur (and take a thickness Dn-i) and of a region where 
n layers occur (and take a thickness Dn) are conceivable 
J76| and have been observed, see e.g. [771 ]. In summary 
of these remarks, we note that it is not uncommon that 
global geometric constraints may destroy the possibility 
of phase coexistence. 

In view of the above discussion, it is of interest also 
in the present case to investigate the use of the (normal) 
pressure p (instead of the strip width D) as the control 
variable. Taking, in the spirit of the general remarks on 
the phase switch method, the appropriate phase switch 
energy cost as an order parameter M, we can sample 
the probability distribution function p{M) which exhibits 
two well separated peaks of generally different weights. 
These peaks are even more clearly visible in the distribu- 
tion of the energy difference p{E^Qrows — E2<jrotus) at fixed 
{u\ as the order parameter M is related to this energy 
difference via a logarithmic function (cf. eq. [T2|) . The 
transition pressure pt is that for which the peaks have 
equal weight (Fig. [9]) and can be determined accurately 
via histogram reweighting. From this we estimate that 
Pt = 22.146 ± 0.015. At the transition, the measured 
misfit A jumps from Ai = 1.913 ± 0.043 (for n = 30) 
to A2 = 1.503 ± 0.046 (for n = 29). Interestingly, the 
misfit where the transition in the NVT ensemble occurs 
(At K, 1.71) is just the average of these two values. 



Comparison of competing candidate stable 
structures 



Returning again to the NVT ensemble, we now con- 
sider the transition from states with 29 layers to states 
with 28 layers. We recall (Fig. 0]) that several different 
candidate structures do exist, and it is not at all clear 
a-prior% which of them should be favored. Again, the 
phase switch Monte Carlo is a convenient tool to solve 
such a problem: we utilize reference states from all four 
of the candidate structures having n = 28 (as shown in 
Fig. |4]) and calculate the free energy difference AF be- 
tween the (unique) structure with n — 29 and these four 
candidates. 

The results (Fig. [T0| clearly show that configurations 
number 1 and number 3 are metastable, because they 
have distinctly higher free energy differences throughout 
the range of A than configurations number 2 and 4 which 
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FIG. 10. Free energy differences between various structures 
with n = 28 rows and the structure with n = 29 plotted vs. 
the misfit A. As configurations nr. 2 and nr. 4 turned out to 
be the same, their free energy curves fall on top of each other. 



practically coincide. In fact, this coincidence between the 
free energies of configurations nr. 2 and 4 is not acciden- 
tal: a closer evaluation of their time evolution shows that 
they transform into each other via sequences of "easy" 
local moves, and although the instantaneous snapshot 
pictures reproduced in Fig. |3] were different, they do not 
belong to different phases in a thermodynamic sense. 

It is also interesting to note that the conclusion that 
structure number 2 is the stable one would not have been 
obtained by a simply comparison of the internal energies 
of the four structures: indeed configuration number 2 has 
the highest energy of all four structures. 

Thus, entropy matters in soft crystals, such as those 
studied here. 



V. CONCLUDING REMARKS 

The principle findings of our study are two-fold: (i) 
We have performed a thorough test of the suitability of 
the phase switch Monte Carlo method for the task of 
determining the relative stability of imperfectly ordered 
structures of typical soft-matter systems, where one must 
deal with systems which have at least one very large lin- 
ear dimension. For such a test, it is crucial to provide 
full information on the model that is studied, and to give 
a careful description of the method and its implementa- 
tion. Moreover we have studied precisely the same model 
system by a thermodynamic integration method thereby 
allowing the first like-for-like comparison between the two 
approaches. We find that the results from both methods 
are compatible, but the accuracy that can be achieved 
using phase switch MC is at least an order of magnitude 
better (Fig. [7]) , despite requiring a factor of 30 less com- 
putational time. 

The reasons for this efficiency gain can be appreciated 
from a glance at Fig. [B] the absolute free energies of 
our system of 3240 particles vary from about 22000 to 



24000 (in suitably scaled units), for a misfit parameter 
A varying from 1 to 2, while the free energy difference 
between the two states that we wish to compare vary 
only from —60 to -|-60 in the same range. These num- 
bers illustrate vividly the basic concept of phase switch 
Monte Carlo: one does better in focusing directly on the 
small free energy difference between the states that one 
wishes to compare, rather than extracting them indi- 
rectly by subtracting two measurements of large abso- 
lute free energies. Thus (in the present context at least) 
phase switch Monte Carlo seems a much more powerful 
approach than thermodynamic integration. In fact, if 
one were to try to bring the errors of the thermodynamic 
integration method down by an order of magnitude - to 
make the error bars of both methods in Fig.[7]comparable 
- one would have to invest a factor of 3000 more com- 
putational time. We feel that the case of relatively small 
free energy differences between competing phases and/or 
structures is rather typical for soft matter systems. In- 
deed for many soft matter systems, such as block copoly- 
mer mesophases, the relative magnitude of free energy 
differences is much less than the factor of about 1/500 
encountered here, and hence such problems could never 
be tackled successfully with thermodynamic integration 
methods since the computational effort to reach the req- 
uisite accuracy would be prohibitive. 

The first problem to which phase switch Monte Carlo 
was applied (in the form of the "Lattice-switch" method), 
evaluated the free energy difference of perfectly ordered 
face-centered cubic and hexagonal close packed crystals. 
Such an application might be regarded as a somewhat 
special case due to the perfect long-range order in these 
defect-free crystals. However, the present work shows 
that the method can equally be applied to imperfectly 
ordered crystals. Here, due to the confinement by struc- 
tured walls together with a misfit between the distance 
between the walls and the appropriate multiple of the dis- 
tance between the lattice rows, somewhat irregular long 
range defect structures form along the walls ("soliton 
staircase"). Additionally several similarly ill-crystallized 
structures can present themselves as candidates for the 
optimal structure (Fig.H]). It would be absolutely impos- 
sible to identify which is the equilibrium structure and 
which structures are only metastable without the phase 
switch Monte Carlo method (Fig. [T0|) . 

We note that the model system that we have chosen 
to study (Fig. [T]) could also be experimentally realized 
in colloidal dispersions, though with some effort: colloids 
coated with polymer brushes experience a short ranged, 
almost hard-sphere-like, repulsive effective potential, and 
bringing them to an interface where water is on top and 
air is below, rather perfect two-dimensional crystals with 
triangular lattice structure form. Interference of strong 
laser fields can be used to create a periodic confining 
potential, through which the misfit and thus the crystal 
structure can be manipulated. We hope that our study 
will solicit some corresponding experimental studies to 
show that the proposed transitions in the number of rows 
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in these crystalline strips actually occur. 

(ii) Our second main finding is that this type of sys- 
tem has an interesting physical property, namely the in- 
equivalence between conjugate ensembles of statistical 
mechanics. When we fix the distance D between the con- 
fining "walls" , the total particle number TV and the total 
(two-dimensional) "volume" V of the system, we realize 
the NVT ensemble. When one studies first order transi- 
tions in the bulk using such an ensemble containing two 
extensive variables (iV, V), a first order transition nor- 
mally shows up as a two-phase coexistence region (e.g., 
at fixed N the two-phase coexistence extends from V/ to 
Vii). However, here such a two-phase coexistence is not 
possible (Fig. [8]) , and thus one has the unusual behaviour 
that at the equilibrium in the "constant Z?" -ensemble the 
conjugate intensive variable (the normal pressure pN^ as 
well as the stress a, cf. Fig. [2]) exhibit jumps (in Fig. [21 
we display the hysteresis loops, but the positions of the 
jumps in equilibrium can be inferred from Ai<" = in 
Figs. [71 and [TUl respectively). When we use a "constant 
p" -ensemble (which is physically reasonable if the con- 
finement of the crystal is effected mechanically in a Sur- 
face Force Apparatus), it is the "volume" (i.e., the dis- 
tance between the walls D) which jumps from Di to Du 
at a well-defined transition pressure, cf. Figs. HI IHl 

One should not confuse this ensemble inequivalence 
with the well-known ensemble inequivalence between 
NVT and NpT ensembles in systems where N is finite: in 
the latter case, the ensemble inequivalence is dominated 
by interfacial contributions (in the NVT-ensemble, when 
Vi < V < Vii, the system is in a two-phase configura- 
tion, as suggested for y — >• oo by the "lever rule", but 
for finite V the relative contribution due to the inter- 
face between the coexisting phases dominate the finite 
size effects). But for T^ — >■ oo these interfacial effects 
become negligible, the properties in the two conjugate 
ensembles are just related by the appropriate Legendre 
transformation. This equivalence between the ensembles 
holds also for liquid-vapor or liquid-liquid unmixing un- 
der confinement in a thin film geometry: when D is finite 
and the particle number N ^ oo, i.e. the lateral linear 
dimensions become macroscopic, we still have ordinary 
two-phase coexistence in the thin films (cf. Fig. [SJ. The 
ensemble inequivalence in the present system arises from 
the lack of commensurability between the thickness D of 
the slit and the appropriate multiple of the lattice dis- 
tance. At a transition pressure pt in the NpT ensemble 
we inevitably have different distances Di, Du between 
the walls for the two phases /, //. Thus, they cannot 
coexist for any uniform value of D. Similar phenomena 
(where the number of layers of a layered lamellar struc- 
ture confined between walls exhibits jump discontinuities 
when D is varied) are already known, both experimen- 
tally and theoretically, for block copolymer mesophases, 
but the aspect of ensemble inequivalence has not been 
addressed, to our knowledge, in these systems studied 
here. 
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VII. APPENDIX 

Here we provide an extended description of the imple- 
mentation of phase switch Monte Carlo, concentrating 
on implementation details at a level suitable for a new 
practitioner. 



ence configuration {M"^} with each phase. The reference 
configuration is an arbitary configuration drawn from the 
set of configurations that are identifiable as 'belonging' to 
phase a. We then associate each particle with a unique 
site of the reference configuration, allowing us to write 



its position r, 
reference site 



^") 



in terms of the displacement Ui from its 



A. IMPLEMENTATION DETAILS FOR THE 
PHASE SWITCH METHOD 

In order to calculate the free energy difference between 
two phases in a single simulation run, the two phases 
have to be linked by a sampling path. In many popular 
approaches, a direct path between the two phases is con- 
structed in the form of a continuous set of macrostates 
associated with the values of some order parameter which 
distinguishes one phase from the other (common exam- 
ples are the total energy or density of a fluid). This path 
traverses mixed phase (interfacial) states Q and is nego- 
tiated using some form of extended sampling to overcome 
the free energy (surface tension) barrier associated with 
the interfacial states. One way to do this is the mul- 
ticanonical method [65| . Alternatively one can directly 
measure free energy differences between successive points 
along the path as is the case in the successive umbrella 
sampling technique |78| . 

In many cases utilizing an inter-phase path that en- 
compasses interfacial states works well, particularly for 
fluid-fluid transitions or lattice models of magnets. How- 
ever, in other cases such a path can be problematic 
Q. For example in the case of solid-liquid coexistence, 
a connecting path will typically run from a crystalline 
phase through several different distinct states including 
droplets of liquid in a crystal, a slab configuration and 
crystalline droplets in a liquid before finally reaching the 
pure liquid phase ^SOj . In such cases the identification of 
a suitable order parameter to guide the system smoothly 
from one pure phase to the other can be difficult, and as 
a result the system may experience kinetic trapping (eg 
in defective crystalline states). 

Thus it is highly desirable to have a method which 
can directly "leap" between the two pure phases (which 
we shall label a, with a = 1,2), avoiding the prob- 
lematic mixed phase states. If the system jumps back 
and forth between these phases a sufficient number of 
times within one simulation run, the relative probability 
with which the system is found in each of them directly 
yields the free energy difference between these phases via 

AF = — In ( p(c=2) ) ■ The phase switch method achieves 

this by supplementing standard local particle displace- 
ment moves (and in the case of a simulation in the NpT 
ensemble, moves which scale the volume of the simula- 
tion box), with moves that switch the system from one 
phase directly into the other phase. This switch is fa- 
cilitated by the representation of particle configurations 
in the two phases. Specifically we associate a fixed refer- 
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Note that whilst there are two reference configurations 
(one for each phase), the phase switch method only con- 
siders one set of displacement vectors which are regarded 
as common to both phases. 

Suppose we are currently in phase a = 1, so that 
the particle coordinates are f^ ' — R- ' + Ui. For lo- 
cal moves in this phase we update particle coordinates 
(in the manner to be described) which, owing to refer- 
ence sites being fixed, is equivalent to updating the dis- 
placement vectors. For a phase switch to phase a = 2, 
we propose a new configuration which is simply formed 
by substituting the reference sites of phase a = 1 with 
those of phase a — 2. Thus the proposed configuration 

(2^ -t ('2) 

is {Fj } — {i?j } + {ui}. If this switch is accepted, i.e. 
if the resulting configuration of phase a = 2 is energet- 
ically acceptable, the simulation will continue to run in 
phase a = 2, again recording the displacements of all of 
the particles from the reference sites of phase a = 2, and 
proposing switches back to phase a = 1 . In this way the 
system switches repeatedly back and forth between the 
phases, allowing one to record the relative probability of 
finding the system in each phase. 

Now generally speaking the displacement vectors that 
characterise phase a = 1 are not typical of phase a = 2 
and thus it will not be energetically acceptable to perform 
the switch from typical configurations of phase a = I. To 
deal with this, one introduces a bias in the accept/reject 
probabilities for local moves that enhances the probabil- 
ity of displacements being generated in phase a = 1 for 
which the phase switch to a = 2 is energetically accept- 
able. The obvious observable to which the bias should be 
administered is a quantity related to the instantaneous 
energy cost of the switch, since this measures how likely 
it is to be accepted. We have employed the switch en- 
ergy order parameter M described in Ref. [22| , which for 
switches from phase a = 1 to a = 2 is defined as follows: 

M(i)^(2)(|^|) = sgn(AE(i)^(2)) . 1j^(^ ^ |Ae(i)^(2)|) 

(12) 
where 



(13) 
where El/f is the energy of the reference configuration in 
phase a, and £^'"'({"1*}) is the energy in phase a, found 
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by applying the displacement vectors {u} to the reference 
configuration {w"^. An obvious substitution gives the 
order parameter for the switch from a = 2 to a = 1. 
Note that an important feature of this definition of this 
order parameter is the logarithm which ensures that the 
binning of the weight function is finer for small values of 
the energy difference and thus serves to ensure that the 
simulation can cover the entire range of M smoothly. 

Now, when implementing local moves for particles, we 
consider not just the energy cost of the move within the 
current phase, but also the change in AI associated with 
the local move via a weight function ri{M). The accep- 
tance criterion for the local move is therefore given by: 

min(l, e-^(E<°)({u'})-E<°)({u}))+r,(M')-,,(M)^^ (.-^^^ 

Note that £;(")({?/})) -£;("({u}) is the energy difference 
due to the move in the phase a that is currently being 
simulated. The energy difference in the other phase is 
only needed for the computation of the new order param- 
eter M' and therefore for the weights r]{M') associated 
with the move. 

Phase switches are generally only accepted from states 
in which M is small - the so called gateway states. One 
instance in which M becomes small is if the displacement 
vectors are themselves small, i.e. if all particles are sitting 
close to their reference positions in both phases. Another 
instance is if there is a high degree of structural similarity 
among the phases, so that the displacements of many of 
the particles in one phase are typical of the displacements 
in the other phase. Note that one does not need to know 
or specify the gateway states to use the method. They 
are sought out automatically when one biases to small 
values of M. 

The acceptance criterion for a phase switch from a — 1 
to a = 2 itself reads: 



a) 



P 



(i)-(2)({4) = niin(l, e-'^(^<''«=»-^<''«"»+-''' '"<'') , 

(15) 
and similarly for the reverse switch. This phase switch 
also includes a weight w to ensure that it occurs with 
a sufficiently high probability in both directions. Note 
that since the phase switch move alters the absolute par- 
ticle coordinates, the associated energy change enters the 
switch acceptance criterion. In the case of phase switch 
simulations in the NpT ensemble, an additional volume 
scaling, must also be taken into account, see below. 

Once suitable weights have been determined (see ap- 
pendix [BJ, one samples the statistics of the two phases 
by accumulating a histogram of the biased order param- 
eter distribution P{M). At the end of the simulation, 
the effects of the weights is unfolded from this distribu- 
tion in the standard manner for extended sampling (9| 
to find the equilibrium distribution P{M). Close to a 
phase transition, this distribution will exhibit two well 



1.5 
0.5- 



1 


- A=1.70 


- 


-- A=1.7r 


- 


i - 
1 ; 


, J 


. , , I 



-20 



-10 



b) 




M 



10 



20 



1.4x10"* 


- 


' 1 ' 1 


- A=1.70 " 


, ^ 4 






-- A-1.71 


g 1.2x10 




1 




o 






,1 


ai.Oxlo"* 


- 


^ 


11 


W 




/ 


11 


' 8.0x10 


- 


/ 


!!• 


1 6.0x10"^ 


- 




h : 


^ 4.0x10 


- 


/ i/i 




o. 

5 




/ '' ' 


i' M\'^\ 


2.0x10 


~ 


/ '' ' 


if v'"' k 






^/ \ , 


, y , ^>^i=.i=:v;" 



-40000 



-20000 



20000 



30 rows 29 rows 



40000 



FIG. A.l. a) The order parameter distribution p{M) for sim- 
ulations at A = 1.70 and A = 1.71 carried out in the NVT 
ensemble, b) For comparison the same distribution is plotted 
against the internal energy difference between the two phases 
for fixed {u\. The order parameter M is deduced from this 
energy difference E2grows — Esorows via the definition given 
in eq. \U\ 



separated peaks, whose areas yield the free energy dif- 
ference as described above. An example is shown in 
Fig. lA.ll (a). Also show in Fig. lA.ll (b) is the distribu- 
tion of the instantaneous energy change under the switch 
E^°' -'({w}) — E'-°'^{{u}) which similarly shows two peaks, 
one for each phase. 

With regard to the choice of reference configuration in 
each phase, in principle this can be an arbitrary configu- 
ration belonging to that phase. In practice, however, for 
crystalline systems one finds that the degree of weighting 
required to access the gateway states can be reduced by 
choosing a reference configuration which is a perfect lat- 
tice. For more general system, eg those with crystalline 
disorder, or for fluids it may be advantageous to try to 
ensure that the particles are not sitting too close to each 
other (eg. by energy minimization of the configuration 
[2l|), since particles which are in close proximity reduce 
the number of gateway states significantly. (We note in 
passing that for fluid systems |22j one requires special 
approaches to guide particles to the gateway states that 
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we will not discuss here as they were not necessary for 
our system.) 



^^6 

15 

f-< 13 
o 
o 
u 

>.12 
11 
10 



4 ' 1 ' 1 ' 1 ' 1 


1 ' 1 ' 1 ' 1 - 


« / -' , 


• ^- new positions . 


_ 


■ ref. positions _ 


-% - ■' A • 


" \ V 


- - ' A \ 


— " % ■v 


: \\ 


1 - % " 


r.^\\. 


\-.-.-; 


n "i ,—| .Ni .'•i 


\, ,%, ,', ,\, - 



104 105 106 107 108 109 110 111 112 
x-coordinate 




-1 1 

displacement in y-direction 



to simply suppress it. Measurements of the distribution 
of displacements in the y-direction is shown in Fig. IA.2I 
(b) and show that preventing particles from fluctuating 
any further in in y-direction than Ay = 0.5 introduces a 
negligible constraint with regard to their natural fluctu- 
ations (and hence on free energy measurements). Doing 
so cured the problem of rare lattice site swaps. 

Finally, we outline briefly how to apply the phase 
switch method in the NpT ensemble. The advantage of 
the NpT ensemble is that the results obtained at one 
pressure can easily be extrapolated to other values of 
the pressure by standard histogram reweighting meth- 
ods. The difference between the NVT and the NpT en- 
semble in this case is that additional volume moves have 
to be carried out in the NpT ensemble. In such moves, 
all particle coordinates are scaled (al ong with the box) in 
both phases in the standard way [y, [22| . Additionally, it 
can prove useful to combine the phase switch move itself 
with a volume scaling move if the equilibrium densities of 
the two phases differ from each other as it was the case 
for our system. For details on the underlying statisti- 
cal dynamics, and acceptance probabilities, see Ref. [2^. 
The problem of particles switching their positions and 
thus creating configurations which prevented any further 
phase switches from being accepted did not occur in the 
case of simulations in the NpT ensemble for our system. 
We obtained (within the error bars) the same free ener- 
gies whether or not we restricted the movement of the 
particles in the y-direction in the way we had to restrict 
them in the NVT ensemble. 



FIG. A. 2. a) A section of a configuration that includes parti- 
cles which have swapped positions in the phase with 29 rows. 
Black squares denote the reference configuration, red (grey) 
lines are the displacement vectors associated with the refer- 
ence positions and red (grey) dots are the new positions. This 
simulation was carried out with 3240 particles at a misfit of 
A = 1.70 in an NVT simulation with 1280 bins for the weight 
function, b) Typical histograms of the particle displacements 
in y-direction. The (almost completely invisible) very small 
peaks at about -1-1 and —1 correspond to particles which have 
swapped their positions. 

With regard to the phase switch simulations of the 
present model of 2d colloids in confinement, we men- 
tion a rare problem that appeared in our simulations of 
the 29 row system. This involved sets of particles on 
neighbouring lattice sites in adjacent rows jumping be- 
tween rows during the simulation, creating in the process 
a ring of particles which occupy each others positions 
(cf. Fig. IA.2I (a)) and remain there. This occurrence 
is primarily a feature of the two-dimensional nature of 
our system, and the well known 'softness' of 2d crys- 
tals. When it occurs it interferes with the operation of 
the phase switch method because the weight function is 
not designed to deal with it, so one is prevented from 
reaching the gateway states. Although one can envis- 
age methods for solving this problem along the lines of 
those used in fluids '221] , our solution to the problem was 



B. IMPLEMENTATION DETAILS FOR THE 
TRANSFER MATRIX METHOD 

The choice of method for determining the weight func- 
tion ri{M) that connects the configurations of high sta- 
tistical weight to the gateway states is to some extent a 
matter of personal taste. A number of approaches ex- 
ist such as the Wang-Landau method [63] or successive 
umbrella sampling [781] . In this work, we have found the 
transition matrix method to be a particularly efficient 
means of determining a suitable weight function. The 
transition matrix method has the advantage that - sim- 
ilar to the Wang-Landau-Sampling - the weights can be 
updated "on the fly" throughout the simulation, allowing 
the simulation to explore an ever wider range of values 
of the order parameter M as the weight function evolves, 
until it eventually encompasses the gateway states of low 
M. Once this has been achieved, one can cease updat- 
ing the weight function and perform a simulation run 
with a constant weight function. An advantage of transi- 
tion matrix method over Wang-Landau sampling is that 
it collects equilibrium data from the outset of the simu- 
lation, whereas Wang-Landau only provides equilibrium 
estimates after a number of preliminary iterations. 

The general idea of the transition matrix method for 
determining weight functions is to record the acceptance 
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probabilities of all attempted transitions and extract the 
ratio of the states' probabilities from it. As all attempted 
transitions contribute to the weight function, including 
those that were rejected, the weight function can be built 
up rather quickly. The details of the implementation are 
as follows and can also be found in [2l|, 122, l6l| and the 
references given therein. 

To implement the transition matrix method, the range 
of the order parameter M, for which a weight function is 
desired, is divided into a number of bins. In our case this 
range corresponds to the values of M that lie between 
the peaks in P{M) which correspond to the two phases 
(cf. Fig. lA.ll (a). A good choice for the binning of the 
order parameter is to choose the bins in such a way, that 
the weight difference between adjacent bins satisfies f22'| 
\ri(Mi^i) — r]{Mi)\ < 2. Then, for every attempted move 
the acceptance probability p, (which is calculated anyway 
for use in the Metropolis criterion) is stored in a collection 
matrix C: 



C{M -^ M') =^ C{M -^ M') +p 



(16) 



At the same time, the probability for rejecting the move 
and thereby keeping the current value of the order pa- 
rameter is also stored: 



C'iM -^ M) => C{M -^ M) + (1 - p) 



(17) 



It is important to note that these probabilities p are the 
"bare" acceptance probabilities and do not include any 
weights. 

The transition probabilities are then simply calculated 
by a normalization of the values in the collection matrix, 
with the sum on the right hand side including all possible 
states to which the system can jump from a given state: 



T{M -^ M') = 



C{M -^ M') 
Efc C{M -^ Mk 



(18) 



In the most general case, this method would create an 
N X N matrix, N being the number of bins or values of 
the order parameter M. In order to derive the correct 
probability distribution from such an N x N transition 
matrix, it is necessary to compute the eigenvector to the 
largest eigenvalue of this matrix. However, it is not neces- 
sarily required to know the exact probability distribution 
in order to create a weight function that will work suf- 
ficiently well. Therefore it is possible to take only those 
transitions occuring between neighbouring bins of the or- 
der parameter into account when computing the weight 
function. In terms of the transition matrix, this means 
that only the diagonal elements - corresponding to tran- 
sitions from a state to itself - and the first off-diagonal 
elements - corresponding to transitions from one state to 
the adjacent ones - are taken into account. Using this ap- 
proach the weight function can be calculated quite easily 
without the need to compute eigenvalues or eigenvectors 
of the transition matrix. In this case, the ratio of the 



probabilities of two adjacent states can be read off di- 
rectly from the transition matrix via 



P{AU 



T{hh 



M, 



i+l) 



yielding the weight difference 

,7(Af,+i)-^(M,) = -ln(^-^^^^+'^ 



(19) 



P{Mi) 



\T{AU 



M^+i) 



Mi) 



(20) 



Of course, when running the simulation, the system is 
still free to perform transitions between any values of M . 
But these transitions are not registered in the transition 
matrix and thus are also not taken into account when 
calculating the weights. In the present study this was 
found to produce accurate and useful weight functions as 
transitions between distant values of M were rare and the 
entries in the second off-diagonal elements of the transi- 
tion matrix were already considerably smaller than the 
ones we used for the calculation of the weights. 
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FIG. B.l. Weight functions for the two-dimensional colloidal 
crystal with A'^ = 3240 particles and structured walls at a mis- 
fit of A = 1.7. The left minimum corresponds to states where 
the system was simulating the phase with 29 rows, the right 
minimum corresponds to 30 rows. Note that the weights have 
an exponential influence on the acceptance criterion. The 
large figure shows the weights plotted against the order pa- 
rameter M as defined in eq. 1121 the inset shows the same 
weight function plotted against the energy difference between 
the two phases in order to illustrate how the definition of 
the order parameter in the logarithm of the energy difference 
stretches the part around M — 0, where phase switches are 
most likely to happen. 

By accumulating the transition matrix in the course 
of a simulation, one obtains an estimate for P{M) which 
can be used to update 'r]{M), thereby allowing the simu- 
lation to explore a wide range of M. Repeated updates 
of 7y(M) thus extend systematically the range of M over 
which one accumulates statistics for the weight function, 
until ultimately one reaches the gateway states. However 
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since updating the weight function during a simulation 
violates detailed balance, we chose to do this at rather 
infrequent intervals of 20000 sweeps. Once the weight 
function extends to the gateway states, we stop updat- 
ing the transition matrix and perform a long phase switch 
simulation with a fixed weight function in order to accu- 
mulate equilibrium free energy data. 

An example of a weight function created for the system 
with TV = 3240 particles (plus 432 fixed wall particles) at 
a misfit of A = 1.7 is given in fig. IB.li also illustrating 
how the definition of the energy order parameter M given 
in eq. [T^l which includes a logarithm of the energy differ- 
ence, leads to a finer binning in the part closer to Af = 0, 
where the phase switches are most likely to happen. In 
fact to ensure that the transition matrix estimate of the 
weight function was sufficiently smooth and reliable in 
this region we reduced the number of bins somewhat. 
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